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The fractional inverse M -7 (real 7 > 0) of a matrix M is expanded in a series of Gegenbauer polynomials. If 
the spectrum of M is confined to an ellipse not including the origin, convergence is exponential, with the same 
rate as for Chebyshev inversion. The approximants can be improved recursively and lead to an iterative solver 
for M 1 x — b in Krylov space. In case of 7 = 1/2, the expansion is in terms of Legendre polynomials, and rigorous 
bounds for the truncation error are derived. 



1. Gegenbauer Polynomials 

The key relation is the generating function 

00 

(\ + t 2 -2tz)^ = Y^t n Cl(z) (1) 

ra=0 

with a real parameter 7 > 0. It defines the Gegen- 
bauer polynomials C%(z), see Q, part 10.9. They 
obey the recursion relation 

{n + l)Cl +1 (z) + {n + 2 7 -l)CP n _ x (z) 

= 2(n + 1 )zCZ{z) (2) 

(n > 0, CZ 1 = 0) and are normalised such that 

T(2 7 + n) 
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Special cases are the Legendre polynomials 

C 1 J 2 (z)^P n (z) (4) 

and the Chebyshev polynomials of the second 
kind 



C\{z) = U n (z). 



(5) 



Estimates of C 7 are obtained with the aid of the 
integral representation 



2 1 - 2 Tr(2 7 + ?i) 

n!r( 7 ) 2 



dip (sin^) 27_1 (z + \Jz 2 - lcos<p) n . (6) 



2. Convergence of Gegenbauer expansions 

The relation (|l|) was introduced to define the 
Gegenbauer polynomials, but it also serves as an 
efficient expansion of (1 + t 2 — 2tz)~~ 1 in polyno- 
mials of z. In this context, the relative error of 
the truncated series is of particular importance: 



R n (z) = 1 - (1 + i 2 -2tzyJ2 tkc K z ) ( 7 ) 
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On the line segment z G [—1,1], the Gegen- 
bauer polynomials are uniformly bounded by 

\CZ{z)\ < C2(l) for ze[-l,l] 

n 2~t-l 

Therefore the expansion (|TJ) converges uniformly 
in z, provided that \t\ < 1. The relative error (||) 
decreases like 

\R n (z)\=OQtr +l ), (10) 

possibly up to powers of n. For the Legendre 
case (7 = 1/2), the following rigorous estimate is 
proven in : 

\Rn(z)\ < \t\ n+1 for ze [-1,1]. (11) 

To extend these considerations to z £ [—1, 1], it 
is convenient to parametrise the complex plane in 
terms of confocal ellipses (as in the case of Cheby- 
shev polynomials, see Pfl). Inside the ellipse 
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z = cosh(6» + i0), 0>O, 0€[O,2tt] 



(12) 



2 



C 7 is bounded by 

\(%{z)\ < C 7 (cosh0) 

< C 7 (l) cosh n9 



(13) 



and convergence of the expansion ([!]) is uniform 
in z if |<|e e < 1. We also have 

\R n (z)\=0((\t\e e ) n+1 ) (14) 

and for 7 = 1/2 the general bound|Q] 

\R n (z)\ < (\t\e e ) n+1 . (15) 

3. Iterative solver 

The aim is to solve 



M^x = b 



(16) 



approximately in Krylov space, i.e. in terms 
of polynomials of M acting on b. To this end, 
parametrise M as 



M = c(l + t 2 - 2tA) 



(17) 



with c, t € C . If M has a spectrum on a line 
away from the origin, we can find a transforma- 
tion with \t\ < 1 such that specA C [—1,1]. More 
generally, if the spectrum of M is bounded by an 
ellipse which does not include the origin, it can be 
mapped into an ellipse (12) with \t\e 6 < 1. In any 
case, |i| resp. \t\e should be as small as possible 
for optimal convergence. 
The formal solution 



x = c~ 7 (l + i 2 - 2tA)' 



(18) 



suggests to insert the Gegenbauer expansion (|l|) 
and to form the approximants 



k=0 



k=0 



The shifts 



s n =c-^CZ(A)b 

inherit the recursion relation 

(n + l)s n+ i + (n + 27 - l)s„_i 
= 2(n + 7 ) As n 



(19) 



(20) 
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(n > 0) to be startet from 

s_i = 
so = c 'b. 

As to the stability of the recursion, the error Ss n 
evolves for large n according to the approximate 
relation 



8s n -i = 2ASs n . 



Consider an eigenvalue A of A and express it as 
A = cosh(^ + i(p) with < 1) < 8. Then the error 
of the corresponding mode of s n behaves like 

Ss n cx e ±{i>+llp)n . 

Note that the shifts s n are to be multiplied by t n 
when accumulated for the solution (|l9|). There- 
fore, the modes of iterated errors in t n 8s n con- 
tribute ~ (|t|e ±1? ) rl < (|i|e e ) n , i.e. error propaga- 
tion is damped due to \t\e e < 1. 

In conclusion, the recursion ( |2l| ) provides a sta- 
ble solver (|l^) for the linear problem (|l^). This 
was confirmed in a first application J5| within the 
boson algorithm for dynamical fcrmions. 

A peculiar feature of this solver is that it cannot 
be started from an arbitrary vector x, but only 
from x-i = or xq = sq as stated above. This is 
due to the fact that M 7 x_i cannot be evaluated 
for arbitrary x-\ (which is needed to bring the 
start vector to the r.h.s. of the equation). If 
an approximant x n is to be improved in further 
iterations, the last shifts s„_i, s n have to be saved 
as well. 

4. Error estimates 

The rest vector and the relative error of x n are 
both controlled by R n eq.(ph: 



b-M^Xn = R n {A)b 
Af-Tr n = R n (A)x 



(22) 
(23) 



Note that, in contrast to the case of plain inver- 
sion (7 = 1), r n is not available during the iter- 
ative process, because M 7 is not computable for 
fractional 7. Therefore mathematical bounds for 
are important to estimate the quality of 
the approximation. (Throughout this paper, || • | 
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designates the Euclidean norm for vectors and the 
induced matrix norm.) 

If Af is normal, i.e. [Af, Aft] = 0, the same 
is true for A, and ||i?„(A)|| is determined by the 
spectrum of A. This results in an estimate for the 
relative error 



\\x-x n \\/\\x\\ < \\R, l+ i(A)\\ 

= max\R n (Xi(A))\. 



(24) 



At this point, uniform bounds like (O) or (15) 
over the spectrum of A are useful. So far, they 
are available for 7 = 1/2 only. 

As an example, let Af = Aft be positive definite 
with specAf c| 
such that eq.([jj 
and Af = \ mm - 



Amm.Amm]- c and t are chosen 
|) maps Af = \ max — ► A = -1 
-> A = 1. This implies 

2 



t 



Jk + 1 



l + t 
1 - 1 



Eq.(10) shows that t is the convergence factor, 
and it agrees with the one for optimal Chebyshev 
inversion and the bound for the Conjugate Gra- 
dient solver. 

If specAf is not known, one may proceed as fol- 
lows: define A, based on a guess about the spec- 
trum, and monitor the norm of the shifts ||s„||. 
The estimate 



\/\\so\\ = ||C2(A)6||/||6|| 

< \\C2(A)\\ 

< C%(coshd) 



(25) 



is assumed to be saturated for large n, and this 
allows for a (pragmatic) determination of 9. If 
\t\e > 1) there will be no convergence and one 
has to stop and try again. Otherwise, is em- 
ployed in the error estimate. It may also help to 
refine the parametrisation of Af (choice of c and 
t) for later use. 



• If the spectrum of Af is known, the rela- 
tive error of the approximation can be es- 
timated. Tight bounds are known for 7 = 
1/2. 

• If the spectrum is uncertain, the growth of 
the shift vectors can be used for a rough 
estimate. 



Polynomial approximations P„(Af) « M" 
are obtained along the same lines. 
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• The polynomials generated by the Gegen- 
bauer expansion are not optimal in any 
sense, but the method is distinguished by its 
conceptional and computational simplicity. 
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5. Conclusions 

• The Gegenbauer expansion allows to con- 
struct iterative solvers for Af 7 a; = b with 
the same rate of convergence as (optimal) 
Chebyshev methods for 7 = 1. 



